Variational Studies of Triangular Heisenberg Antiferromagnet in Magnetic Field 
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We present a variational study of the Heisenberg antiferromagnet on the spatially anisotropic 
triangular lattice in magnetic field. First we construct a simple yet accurate wavefunction for the 
1/3- magnetization plateau uud phase on the isotropic lattice. Beginning with this state, we obtain 
natural extensions to nearby commensurate coplanar phases on either side of the plateau. The 
latter occur also for low lattice anisotropy, while the uud state extends to much larger anisotropy. 
Far away from the 1/3 plateau and for significant anisotropy, incommensurate states have better 
energetics, and we address competition between coplanar and non-coplanar states in the high field 
regime. For very strong anisotropy, our study is dominated by quasi-ld physics. The variational 
study is supplemented by exact diagonalization calculations which provide a reference for testing 
the energetics of our trial wavefunctions as well as helping to identify candidate phases. 
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I. INTRODUCTION 

The spin-1/2 Heisenberg antiferromagnet on a spa- 
tially anisotropic 2d triangular lattice is a deceptively 
simple spin system which nevertheless possesses very 
rich physics. Despite having attracted much attention, 
a complete understanding of the model has not been 
achieved. This can be attributed to the enhanced quan- 
tum fluctuations arising from a combination of low di- 
mensionality, small spin, geometrical frustration and spa- 
tial anisotropy, thus leading to a rich phase diagram. At 
zero field, studies have suggested that the anisotropic 
system may possibly remain disordered even at zero tem- 
perature. In experimental realizations of the triangular 
antiferromagnet, a 1/3-magnetization plateau was found 
for the approximately isotropic material Cs2CuBr4 but 
not for the more anisotropic Cs2CuCl4ii~— A recent ex- 
perimental study of Cs2CuBi'4 further revealed a cascade 
of phases in the fields above the 1/3 plateau, which are 
still not understood, 7 

Analytical studies on the model have been done 
on specific regions of the phase diagram, for instance 
low anisotropy near the 1/3-magnetization plateau*^ 
large anisotropy limit ^ and high field limit iliiiS Sev- 
eral numerical studies using exact diagonalization^ - — 
series expansion) 17 ' 18 coupled-cluster method^ den- 
sity matrix renormalization group ; 20 ' 21 and variational 
approaches^— have also been used to analyze the 
model. Motivated by the experimental and theoretical 
works, we perform a variational study using simple yet 
powerful wavefunctions, attempting to cover a large por- 
tion of the phase diagram. 

We consider the Heisenberg model in external mag- 
netic field h, 
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where S r is the spin operator on site r and J rr i are the 
nearest neighbour exchange couplings. Throughout, wc 




FIG. 1: Triangular antiferromagnet with coupling constants 
J and J' between nearest neighbours along horizontal and 
oblique directions respectively. Three sublattices A, B, and 
C important in the isotropic case are also labelled. 
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FIG. 2: Spin orderings on the isotropic triangular lattice in 
the field, where three arrows refer to three sublattices in- 
dicated in Fig. [1] a) spiral (non-coplanar umbrella); b,c,d) 
coplanar Y, uud, and V. 



extensively use hard-core boson picture, mapping = 
b r and S* — h — n r : 

H = — t rr i (b\b r i + h.c.) + J rr in r n r i (2) 

{rr') {rr') 

—jj, n r + const , (3) 



1 u 1 ST t 

t rr i — — J rr fj [1 — 2 ^ Jrr' 

j.1 (Z r 
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The boson hopping amplitudes are negative and there- 
fore frustrated on the triangular lattice, making this a 
challenging interacting problem. 

Figure [2] depicts spin ordered phases considered in our 




FIG. 3: Boson interpretation of the spin orderings in Fig. [2] 
Gray scale shows density order (n r ), while arrows show su- 
perfluid order (&£). a) Spiral is a uniform superfluid phase 
with rotating phase angles; c) uud is a Mott insulator with 
V3XV3 CDW order; b,d) Y and V are supersolid phases 
which contain both charge density and particular superfluid 
orderings. 



variational study. While these are simple to draw, re- 
alizing them as wavefunctions is non-trivial. The spiral 
phase is a boson superfluid containing rotating phase an- 
gles (see Fig. [3^) and is captured by an elegant Husc-Elscr 
generalization of the Bijl- Jastrow wavefunction.— For the 
other phases, constructing simple and yet accurate wave- 
functions is not as straightforward and requires consider- 
ation of their physical nature in terms of bosons. Thus, 
the uud is a Mott insulating phase (see Fig.[3fc) which re- 
quires a wavefunction with strongly localized bosons and 
rapidly decaying correlations. The Y phase is an interest- 
ing supersolid phase (see Fig. [5p) with rapidly decaying 
boson correlations between sites on one of the three sub- 
lattices as well as long-range correlation between sites on 
the other two sublattices. The V phase is a different su- 
persolid (see Fig. with long-range boson correlations 
between all sites; here we find that two different construc- 
tions of the trial wavefunctions are required to capture 
the lower and higher density regimes. (We remark that 
supersolid phases of bosons on the half-filled triangular 
lattice and in the presence of stro ng rep ulsion have been 
of much recent interest, see Refs. 

and citations 

therein.) 

In section [Til we present simple, few-parameter can- 
didate wavefunctions used in our isotropic study. En- 
couraged by the accuracy of these candidates, we gener- 
alize the wavefunctions to incommensurate versions for 
our anisotropic study in section UTTl We conclude with a 
discussion of the results and implications for Cs2CuBr4 
in section II VI 
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TABLE I: Comparison of uud trial energies for different num- 
ber of short-range Jastrow parameters for 12 bosons on a 6 x 6 
cluster. Localized orbitals in the permanent extend only to 
nearest neighbour sites with amplitude a, which is single vari- 
ational parameter in the first listed uud case. The second uud 
case has one nearest-neighbor (nnb) Jastrow pseudopotential 
which is taken to be the same between any pair of nnb sites, 
while the third case has two such parameters, one for A — B 
and A — C nnb pairs and the other for B — C nnb pairs, as 
is appropriate given lattice symmetries of the uud state. We 
also show trial energy for the spiral state with 4 variational 
parameters (same as in table ITTT|) ; this state performs poorly 
compared to the uud state. 



II. ISOTROPIC TRIANGULAR 
ANTIFERROMAGNET: 6x6 STUDY 

In this section, we consider the isotropic triangular 
Heisenberg antiferromagnet with J rr i = J for all nearest 
neighbour links. Beginning at density n = 1/3, where 
we have an excellent wavefunction for the uud Mott in- 
sulator phase of Fig. [3J:, we construct similarly inspired 
wavefunctions for nearby supersolid phases of Figs. |3j3,d. 
Next, we describe Bijl-Jastrow-type wavefunctions for 
the spiral of Fig. [3^ and V supersolid of Fig. |3ji, which 
perform better for n further from 1/3. We also discuss an 
alternative construction of the uud state using a det x dot 
("2-parton") trial wavefunction. For the ED calculations, 
we compute ground state energies for Nb < 12, while 
Ni, > 13 data is taken from Bernu et al^ 

We perform all studies at fixed boson number Nb. For 
each wavefunction below, we also include a Jastrow factor 

Jastrow({n r }) = e"^r/ u rr ,n r n r , (5) 

with simple choices of pseudopotentials u rr i providing 
additional variational freedom. 



A. uud state at n= 1/3 

At density n = 1/3, the uud phase is stabilized 
by quantum fluctuations. We construct a simple bo- 
son wavefunction by using Nb = N/3 orbitals localized 
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TABLE II: Comparison of Y trial energies for different num- 
ber of short-range Jastrow parameters for 13 bosons on a 6 x 6 
cluster. The Y state is constructed by adding one boson to 
the uud state as described in the text; we allowed the same 
variational parameters as in the uud case in Table [I] We 
also show trial energy for the spiral state with 4 variational 
parameters, which has higher energy than the Y state. 



around sites Aj, j = 1 . . . N/3, from sublattice A: 



IVw) 



<# c (r) 



({r fc }|Vw) = Perm [^(r,)] 



N/3 / \ 

n (e^ oc m 6 m i°) > ( 6 ) 

1, r = A K 

—a, r = neighbour of Aj , (7) 
0, otherwise 

(8) 



For a = this reduces to the classical CDW state with 
bosons strictly localized on the A sublattice and mini- 
mizing the potential energy. Non-zero a allows bosons 
to hop to nearest-neighbor sites and gain some kinetic 
energy; a > is appropriate for boson hopping t rr i < 0. 
In Eq. ([5]), column j of the Permanent matrix is given by 
the j-th orbital (centered on Aj) evaluated on the occu- 
pied sites |ri} i 27 ' 28 One can loosely connect this wave- 
function with a picture starting from the "Ising" limit, 
J z 3> \t\, and pcrturbatively building in boson kinetic 
energy effects.— 

Table [J compares the trial energies for different num- 
ber of variational parameters. Excluding any Jastrow 
factor, the single-parameter trial state already captures 
the important exchange energies; for example, it is closer 
to the exact ground state in the zero momentum sector 
than to the first excited state in this sector (not shown) J£ 
Adding a short range Jastrow factor further improves the 
trial energy, while longer range Jastrow parameters are 
unimportant since correlations in the Mott insulator de- 
cay rapidly (see Fig. [11] in Appendix IHj). We see that 
the simplest localized orbitals extending only to near- 
est neighbour sites (and with relatively small amplitude 
a ~ 0.23) perform very well, which suggests strong uud 
order in the 1/3-fillcd systenii&i 7 - Indeed, for the optimal 
wavefunction, we calculate the boson density to be 0.76 
on A sites and 0.12 on BC sites. 



we construct a candidate for the Y supcrsolid phase by 
adding bosons to an extended orbital on BC: 

I v N b ~N/3 

\^y) = (X>fc(r)4j IVw) , (9) 

-1, r e B 

<t>Tc{r) = { -1, r e C (10) 
v 0, re A 

Just as in the uud case, the wavefunction can be written 
as a Ni, x Nf, permanent. The first A^/3 columns con- 
tain the same <fi l ° c orbitals as in the uud state, while the 
remaining Nb — A^/3 columns all contain the extended 
orbital residing on the B and C sublattices. The 
alternating signs of the extended orbital are appropriate 
for bosons hopping on the BC honeycomb with t rr > < 0. 
Nearest neighbour contacts on BC are suppressed by 
adding a Jastrow factor. 

Table [II] compares the Y energy against spiral and ED 
energies for 13 bosons on the 6x6 cluster. Our Y state 
is close to the ED ground state from Bernu et al^ thus, 
the trial energy is below the first excited state in the 
same sector (not shown), while the spiral is significantly 
higher. 

We consider such Y states for all boson densities above 
1/3 and find them to give lowest trial energies among all 
our states for Nt, = 13, . . . , 15. We discuss properties of 
the Y states in Appendix [B] Here we note an interest- 
ing feature that boson correlations are long-ranged for B 
and C sublattice sites but are short-ranged for A sites. 
The A sublattice remains "Mott-insulating" despite the 
superfluid on the BC honeycomb. The absence of the 
"proximity effect" on the A sublattice is due to cancella- 
tions from alternating superfluid order parameter on the 
B and C sublattices^ In particular, just as in the uud 
case, we cannot construct Bijl- Jastrow-type wavefunction 
for the F-state. 



C. Spiral state at n < 1/2 

At half-filling, the 120° magnetically ordered state 
(spiral) is believed to be the ground state. We use Huse 
and Elser wavefunction^ which generalizes Bijl-Jastrow- 
type wavefunction by including complex 3-body terms, 
to accurately describe the corresponding superfluid state 
of bosons near half-filling. In this wavefunction, all the 
bosons reside on an extended orbital, 



N h 



• ? bt 



|0>, (11) 



B. Y state at n > 1/3 

Starting from the uud wavefunction where we have 
2;ood exchange energies between sublattices A and BC, 



with Q = (An/3,0). Despite frustration, the bosons gain 
some kinetic energy while hopping along any lattice link. 
Nearest neighbour contacts are suppressed by adding a 
long range Jastrow factor. The three-body phase fac- 
tor, which respects the symmetries of the classical state, 
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TABLE III: Comparison of spiral energies for different num- 
ber of variational parameters for 18 bosons on a 36-site clus- 
ter. The extended orbital of the spiral state has amplitudes 1, 
e 12 ^ 3 , and e l47r//3 for sites on sublattices A, B, and C respec- 
tively. The first spiral case has one nnb Jastrow pseudopo- 
tential which is taken to be the same between any pair of nnb 
sites, while the second spiral case has an additional parame- 
ter 7 for the Huse-Elser phase factor. The third case has two 
more parameters w and p for a long range pseudopotential 
|iJ^| P between any pair of sites i and j. We also show trial 
energy for the Y state with 4 variational parameters, which 
has slightly higher energy than the spiral state. 



serves as an additional variational parameter. For de- 
tails, the reader is referred to the original Ref. [22|. 

Among our trial states, the Huse-Elser wavefunction 
has lower energy than the Y state for Nb = 16, ... , 18, 
but only by a very small amount (cf. Fig. [4]). Table [TTT1 
shows that the 18-boson spiral energy is only slightly 
lower than the Y energy. This is perhaps not surpris- 
ing since the classical 120° order may be viewed as the 
spiral or K-shape order depending on the plane's orien- 
tation. A recent variational study^i using different con- 
structions of the spiral and Y states obtained —0.1827 
for their many-parameter spiral state, which is also lower 
than their Y trial energy by a small amount similar to 
that in our study. Other recent work o 21 i 26 observed an 
abrupt change from the spiral to Y supcrsolid in the half- 
filled model as the spin anisotropy is varied through the 
SU(2)-invariant Hciscnbcrg point. In principle, thinking 
in terms of wavef unctions, the spiral and Y can be dis- 
tinct phases with different postulated symmetry breaking 
also in the SU(2)-invariant model. However, this could 
also be a plane reorientation transition, and the closeness 
in energy of the Y trial states reflecting their ability to 
capture the 120° spiral order. 



D. 



V state at n < 1/3 



Let us now consider densities slightly less than 1/3. 
We start from the uud state and picture it as a filled A 
sublattice. An appealing scenario is to introduce holes 
and let them move around on A and condense. We auto- 
matically retain charge order selecting the A sublattice 
vs B and C. The condensation of holes on the A gives 
boson supcrfluid order there and by proximity effect also 
on the B and C sublattices. Since tAB = t-AC < 0, we 
expect the phase angle on the BC to be shifted by it from 
the A. The resulting supersolid is precisely the V state. 

Direct wavefunction implementation of this scenario is 



TABLE IV: Comparison of V^erm energies for different num- 
ber of variational parameters for 11 bosons on a 36-site clus- 
ter. The Vpcrm state is constructed by removing one boson 
from the uud state as described in the text, and the param- 
eters here are of the same type as in the uud case in Table 
I. We also show trial energy for the spiral state with 4 vari- 
ational parameters, which has higher energy than the Vperm 
state. 



described in Appendix [X] and leads to a sum of perma- 
nents, which becomes prohibitively costly to evaluate for 
more than few holes. In the appendix, we also motivate 
a qualitatively similar wavefunction with a simpler am- 
plitude given by a single permanent, 



{{rk} |Vv p , 



= Perm 



/ #>°(n) 



P° c (r Nb 



V 1 



1 



1 / 



.(12) 



The first Nb rows contain the localized orbitals of the 
uud construction evaluated at the boson positions, while 
the remaining N/3 — Nb rows are filled with 1-s corre- 
sponding to "zero wavevector" condensate of holes (see 
Appendix [S] for details) . 

Table |IV] shows the V pe rm energy for 11 bosons on 
the 6x6 cluster. Being a descendant of the excellent 
uud state, even with no Jastrow factor the Vperm per- 
forms very well and lies roughly half-way between the 
ground state and the first excited state with the same 
quantum numbers (the latter is not listed in the ta- 
ble). In particular, the V perra clearly wins over the spi- 
ral superfluid with uniform density. Just as in the uud 
case, adding simple short-range Jastrow parameters fur- 
ther improves the trial energy of the Vp er m state. At 
this stage, we did not include long-range Jastrow pseu- 
dopotcntials, which would be needed for correct long- 
wavelength description 2 ^ of superfluid correlations in the 
V phase. 

The Vp erm state gives our best variational energies for 
N b = 6, . . . , 11. In Appendix E we measure properties 
of this state and verify the supcrfluid order with opposite 
signs on the A and BC sublattices as anticipated above. 
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TABLE V: Comparison of Vbijl 
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trial enerj 
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ferent number of variational parameters for 3 bosons on a 
36-site cluster. The extended orbital of the VBiji- jastrow state 
has amplitude e M ^ 2 on sublattice A and —1 on sublattices B 
and C; fi is the only variational parameter in the first listed 
VBiji-jastrow case. The second VBiji-jastrow case has one nnb 
Jastrow pseudopotential which is taken to be the same be- 
tween any pair of nnb sites, while the third case has two ad- 
ditional parameters for a long range pseudopotential which 
is the same as in the spiral case in Table IIIII We also show 
trial energy for the spiral state with 4 variational parameters, 
which has slightly higher energy than the Vayi- jastrow state. 

E. V state at n < 1/3 

The above wavefunction for the V phase is obtained 
from the strong uud state and a priori is not expected to 
remain good at low density. Here we consider an alterna- 
tive construction of the V supersolid using Bijl-Jastrow- 
type wavefunction, working directly with bosons and con- 
densing them into an appropriate extended orbital, 

I VV.Bijl- Jastrow) = ^E^W 6 ^ I ) ( 13 ) 

^(r) = { ^V^C ^ 

This orbital has opposite signs on the A and BC sublat- 
tices as expected from Fig. [3ji. The "chemical potential" 
yit on the A sublattice allows us to control the charge or- 
der. Similar to other wavefunctions for states with super- 
fluid order, it is necessary to include a long-range Jastrow 
pseudopotential. 

As we argue below, this wavefunction is a natural can- 
didate at low boson densities. On the 6x6 cluster, it 
optimizes better than the V pcim state for Nb < 6 and 
also has better energy than the spiral state, see Fig. [4] 
As an example of variational results, Table IVl shows the 
VBiji- jastrow energy for 3 bosons on the 6x6 cluster. We 
see that the V state is slightly better than the spiral state. 
However, both states are quite close in energy and close 
to the exact ground state. We discuss this more below 
and see what we can infer about the competition between 
the coplanar and spiral states from ED spectroscopy. 

First, we want to connect the competing Psiji-jastrow 
and spiral states with physics at low boson densities. 
In the absence of interaction, the kinetic energy mini- 
mizes at two distinct points in the Brillouin zone, ±Q = 
±(47r/3, 0). Boson condensation at one point gives rise 
to the spiral phase; schematically, the spiral wavefunc- 



tion is given by (b-^lO) [or (&' ^) Arb |0) for the oppo- 

site wavevector]. A more complex condensation pattern 
including both points produces a coplanar state, with 
schematic wavefunction (e M 6^.+h.c.) Arb |0). When a = 0, 

this gives boson orbital </>(?•) = cos(Q ■ r) taking values 
{+1,-1/2,-1/2} on the three sublattices, which is es- 
sentially the (fry* orbital in Eq. (fTfj) . On the other hand, 
a = 7r/2 corresponds to a different state with zero bo- 
son density on one sublattice and alternating supcrfluid 
phases on the other two sublattices; in terms of spins, 
this is a coplanar "'J" -type state which has similar sym- 
metry to the Y state in Fig[2}D, but with the vertical spin 
flipped up. For either V or VP, there are two more degen- 
erate states given by lattice translations or equivalently 
by adding ±27r/3 to the phase a. ReferencelTTI studies the 
dilute boson problem analytically for the isotropic lattice 
and predicts that four-boson interactions select coplanar 
states; this study does not resolve between the V and VP 
states, which requires considering six-boson terms. 

Returning to our example with 3 bosons on the 6x6 
lattice, the ED ground state has momentum quantum 
number k = Q (there is also a degenerate state with 
opposite momentum), while we also find two very close 
states with k = 0, one with even and the other with odd 
parity under inversion. This can be traced to 4 degener- 
ate eigenstates of the kinetic energy, 

Our spiral wavefunction construction gives essentially the 
first two states with momentum quantum number k = 0. 
Our 3 degenerate V-states, upon making combinations 
that are momentum eigenstates, give an even-parity k = 
combination as well as the last two states with k = ±Q 
from Eq. (fl"5j) . Finally, the 3 degenerate VP-type states 
give an odd-parity k = combination and the same two 
k = ±Q states. It is clear that these trial states are not 
independent for this small number of bosons; we cannot 
resolve the phases, but we can start looking for some 
tendencies. For example, we can view the fact that the 
k = ±Q are lower in energy than k = as an indication 
for the coplanar states being better than the spiral. In 
principle, we could also try to resolve between the V and 
^ by comparing the even/odd-parity k = states, but 
the splitting is too tiny. 

We have similarly examined ED spectra with Nb = 
4, . . . , 9 bosons, paying attention to near degeneracy of 
ground states and their quantum numbers. The resolu- 
tion between the spiral and coplanar states due to inter- 
actions becomes clearer with increasing density, and in 
each instance the ED data is consistent with the copla- 
nar states being better, which is also supported by the 
VMC data. As far as the resolution between the V and 
'J states is concerned, we can not tell anything with bo- 
son number below 6, while for higher boson density we 
start seeing evidence in favor of the V state. The V 
state is expected coming from the n = 1/3 plateau as 
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FIG. 4: Comparison of variational energies (per bond) for the 
36-site cluster with isotropic exchanges. The trial energies for 
the planar states are obtained using Bijl-Jastrow-type wave- 
function for V (Nb < 6) and permanent-type wavefunctions 
for V (7 < N b < 11), uud (N b = 12), and Y (N b > 12). The 
spiral trial energies are obtained using the Huse-Elser wave- 
function. For clarity, ED ground state energies are subtracted 
at respective boson numbers. The classical energy curve pro- 
vides a reference for judging stabilization of specific phases 
by quantum fluctuations. 



we discussed earlier. One possibility is that V occurs for 
all n < 1/3, but we cannot rule out transition to the ^/ 
coplanar state at low densities. Our VMC study of the 
spatially anisotropic model on larger clusters in Sec. lIII El 
also suggests that the coplanar phase (incommensurate 
in this case, so there is no distinction between V and 
\E') wins over the spiral also for a range of anisotropics, 
strengthening the conclusions here on the coplanar versus 
spiral energetics. 



v 




1 2 3 4 5 

Magnetic Field h/J 

FIG. 5: Magnetization curve of the 36-site cluster obtained 
using the variational energies. The pronounced uud plateau 
agrees well with the ED results in Ref. [Til . This shows that 
the variational tool is able to capture this Mott-insulator as 
well as nearby supersolid phases. 



G. Magnetization process on the isotropic lattice 

Using the trial energies from our studies at fixed Nb, 
we can work out the magnetization curve as a function 
of field h. Fig. [5] shows this for the 36-site cluster. The 
boundaries of the 1/3 magnetization plateau are deter- 
mined by the energy gaps to adding or removing one 
boson to the uud state. Since our permanent construc- 
tions give very good trial Y and V states in this regime, 
the estimate of the plateau range is quite accurate. To 
check finite size effects, we repeated the calculation on 
a 63-sitc cluster and obtained critical fields H c i w 1.4 J 
and H c2 « 2.2J. 



H. 2-parton trial wavefunctions and alternative 
construction of uud state at n = 1/3 



F. Summary of trial energies on the isotropic 
lattice 

Fig. [4] summarizes the spin energies (per bond) of com- 
peting trial states calculated for the 6x6 cluster with pe- 
riodic boundary conditions for all A^. For a better com- 
parison of the accuracies of these trial wavefunctions, we 
subtract the ED ground state energy at each boson den- 
sity. The energies of the classical state are included to 
emphasize the stabilization of specific phases by quantum 
fluctuations. Our wavefunctions are particularly accurate 
in the vicinity of the plateau, and also at low boson den- 
sities (higher fields). In the latter regime, the classical 
energies approach ED values at low densities, indicating 
vanishing quantum fluctuations^ ' 11 ' 12 



The above direct study using spiral, Y, uud, and V 
states is sufficient to describe the phase diagram of the 
spatially isotropic triangular antiferromagnet in the field. 
We now consider a versatile set of trial wavefunctions 
which we will call "2-parton" states. One motivation is 
to give a practical realization of the Chern-Simons flux 
attachment treatment in Ref. l30l . (The relation between 
the par ton and Chern-Simons approaches is discussed in 
Ref. l3l| and citations therein.) Another motivation is to 
prepare for an anisotropic lattice study in Sec. IIIII We 
should say from the outset that while such parton con- 
struction is typically used to produce fractionalized (spin 
liquid) states, it can also be used to give more conven- 
tional states such as CDW of bosons with no topological 
order as discussed below. 

We represent the boson operator in terms of two 
fermions, b = d^S 2 \ subject to constraint b^b = 
c ^( 1 )t c ;(i) — d( 2 )tjf( 2 ) n each site. Imagine some "mean 
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field Hamiltonian" for each parton flavor, 

= -E(i^i^'^ )t 4" ) +^-) 

(rr') 
r 

Here we write the parton hopping amplitudes (which can 
be complex) as = |t^|e*°rr'; we also allow site- 

(n) 

dependent chemical potentials rir to test CDW tenden- 
cies. The 4"' and 

/ir™' are variational parameters. We 
solve and fill up the corresponding Fermi seas with 
Nd-t = Nd 2 = Nb particles. A valid bosonic wavefunction 
is obtained by applying a Gutzwiller projection such that 
every site is either empty (rib = nd 1 = rid 2 = 0) or con- 
tains both g?i and di partons (rib = ll d t = n<i 2 = 1) : 

\^2 P ) = Pg n 4x )+ 4f + i°>- ( 17 ) 

q n <£FS„ 

For each boson configuration, the amplitude is given by 
a product of two Slater determinants. One feature of this 
construction follows from the fermionic statistics which 
provides an inherent repulsive Jastrow effect for the par- 
ticles. This effect can be tuned as follows, 

({r k } \^2 P ) = deti ■ det 2 ■ Idet^ 1 " 1 • Ide^p" 1 , (18) 

which preserves the "sign structure" of the wavefunction 
while allowing more variational freedom with parameters 
Pi and P2 ■ Numerical calculations can be performed us- 
ing efficient determinantal Monte Carlo techniques^ 

Besides treating boson repulsion, we want to have good 
kinetic energy. We can write the frustrated boson hop- 
ping amplitudes in Eq. <(4]) as z®, — \t^,\e za r-r' and view 
this as a problem in an external orbital field producing 
flux 7r through each triangle^ To capture this in the par- 
ton treatment, we view d^ and d^ as charged particles 
whose charges add up to that of the boson b; we therefore 
require 



Thus the parton mean field Hamiltonian should contain 
fluxes such that for the two flavors they add up to the 
original flux seen by the bosons. We can still make dif- 
ferent choices, say, for the di; however, once the a^J, are 

(2) 

fixed, then the uniquely determined. 

We first discuss what we will call "Chern-Simons" 
states that realize the idea in Refill For the d\ hopping, 
we take uniform flux of nir per triangle, where n is the bo- 
son density per site. With this choice, the deti Slater de- 
terminant fills the "lowest Landau level" band and gives 
a finite lattice version of the usual Chern-Simons factor 
rii<j ( z i~~ z j)- We can loosely view the det! as performing 
flux attachment transformation from the bosons to the 
di fermions^i Upon subsequent "flux smearing" mean 



field, the di see flux (1 — n)ir per triangle. In the ab- 
sence of site-dependent chemical potentials and for some 
rational densities, the det 2 Slater determinant is gapped, 
and the boson wavefunction realizes a fractionalized "chi- 
ral spin liquid" j 30 i 33 ' 34 We have tried these "topological" 
states for several densities such as n = 1/3, 1/4, 1/6 on 
the isotropic triangular lattice and found that they are 
poor compared with the uud and V states described ear- 
lier. Thus the interesting proposal of plateaus due to 
chiral spin liquid states is not realized on this lattice.— 

We now specialize to density n = 1/3 and allow a 
chemical potential on the A sublatticc: [Ia ^ 0,/is = 
He = 0. We find that optimal fj/^ , ii^ are large and pro- 
duce strong CDW order in the mean field state. When 
this happens, the trial boson state Eq. (fT8"]) is no longer 
topological in nature. Indeed, if the parton hopping is 
set to zero, this construction simply gives the classical 
V3 x Vo CDW state. The particles completely occupy 
sublattice A, and there is a large gap at the parton Fermi 
level. Adding small hopping does not close this gap but 
only builds in some charge fluctuations into the parton 
mean field and thus into the boson trial state. Working 
perturbatively in / , the leading modification to 
the classical boson CDW wavefunction is to add config- 
urations where one particle moves from a site Aj to a 
neighbor r. The amplitude for such a configuration is 

proportional to ^.r^r/O^A ) ~ ^A^r/^A ' ' wnere 
we have kept track of all signs and introduced schemati- 
cally boson charge gap /x^ . The result is similar to the 
perturbative picture of the CDW working directly in the 
boson language that motivated the wavefunction Eq. (|6|) . 
Thus at this level the 2-parton states with strong CDW 
potential are qualitatively the same as the permanent 
uud state in Sec. Ill Al 

The above leading order structure holds for all 2- 
parton states satisfying Eq. (fTT)]). At higher order, the 
states will differ, and amplitudes can be complex in 
general: e.g., the Chern-Simons wavefunction described 
above is complex- valued. On the other hand, the perma- 
nent uud wavefunction is real. The boson Hamiltonian 
Eq. @ is invariant under complex conjugation in the 
number basis, and the uud state preserves this symme- 
try. We can construct a real-valued 2-parton state by 
taking the d\ fluxes to be or 7r through up or down 
triangles; the di partons see correspondingly it and 
fluxes. We will call this state U U1B" . It was origi- 
nally discussed at half-filling in Ref. HH, where (in the 
absence of chemical potentials) it has Dirac nodes at the 
Fermi level and realizes so-called Algebraic Spin Liquid 
state. This particular state has a good trial energy in the 
Heisenberg model ) 24 ' 25 i 31 and can be viewed as a more 
elaborate real-valued version of the Laughlin-Kalmeyer 
state (see Sec. IIC of Ref. [3l| for more discussion). Away 
from half-filling, the U1B mean field state has Fermi sur- 
faces of partons and may be unstable to a mechanism 
described in Ref. l36l However, this is not a direct con- 
cern here since we are gapping out the state by adding 



8 



large \xa potential and are connecting to the strong CDW 
of bosons. The virtue of using the 2-parton framework 
is that it naturally builds in small charge fluctuations 
as described above, and determinants are easier to com- 
pute as opposed to permanents. Using this construction 
for the isotropic 6x6 lattice at n = 1/3, we obtain a 
very competitive energy —0.1341 (cf. trial energies in Ta- 
ble IT]) with [la ~ 2 and p sa 0.75. (We also obtain close 
trial energy using the Chcrn-Simons state with strong 
CDW potential, in agreement with the earlier discussion 
that all 2-parton states can similarly capture leading lo- 
cal charge fluctuations when the charge order is strong). 
The 2-parton constructions are particularly useful on the 
anisotropic lattice to be discussed in Sec. IIII| since they 
naturally connect to the decoupled chains limit and al- 
low us to detect where the quasi-ld physics sets in and 
explore CDW instabilities. 



III. ANISOTROPIC TRIANGULAR 
ANTIFERROMAGNET 

Motivated by the unknown phases of Cs2CuBr4 in the 
fields we extend our study to the spatially anisotropic 
lattice. As the phase diagram is much more complex, it 
is appropriate to begin this section with a short review of 
the different regimes and phases discussed in theoretical 
literature. Next, we describe anisotropic extensions of 
the wavefunctions introduced in Sec. Inland then present 
our variational results. 



A. Review of phases from theoretical studies 

In this review, it is convenient to refer to a schematic 
phase diagram shown in Fig. 1101 where we parametrize 
the anisotropy using 5=1 — J'/ J. 

In zero magnetic field (bottom axis in Fig. I10[) . varia- 
tional studies suggest that the spiral phase remains stable 
for small lattice anisotropy, while the more anisotropic 
region may contain one or two spin liquid phases i 22 ' 23 i 25 
This is supported by an ED/DMRG study which found 
signatures of spin liquid for J' /J < 0.78 from numer- 
ical measurements of spin structure factor, excitation 
energy gap, and spin correlation^ However, for large 
anisotropy, an analytical study near the decoupled chains 
limit predicts a collinear antiferromagnetic order This 
suggests that the zero field limit is a challenging region 
which remains unsettled. 

In the high field limit near full polarization (top phase 
boundary in Fig. [TD|) . analytical studies of the dilute bo- 
son gas show that the V phase (commensurate and in- 
commensurate) is the likely candidate near the spatially 
isotropic regime jii while the spiral phase dominates for 
strong anisotropy.— 

At intermediate fields, an interacting spin wave ex- 
pansion about the uud plateau in large S and low 
anisotropy limit (left axis in Fig. ITUl) shows that the 



plateau extends considerably into the anisotropic region, 
with commensurate coplanar Y and V phases present 
next to the plateau— In addition, incommensurate copla- 
nar and distorted spiral phases are also predicted for 
larger anisotropy. 

In the nearly decoupled chains limit (right axis in 
Fig. HQ), Ref-Hfl argues that the interchain coupling is a 
relevant perturbation and can induce various boson CDW 
phases or a spiral phase. The former happens for small 
and intermediate fields, while the latter is expected near 
the saturation field. 



B. Anisotropic versions of wavefunctions 

Although our study began with a goal of construct- 
ing accurate wavefunctions for identifying the unknown 
phases of Cs 2 CuBr4, it quickly became clear that this 
is far from easy. The number of theoretically possible 
phases reviewed above is already very rich, with different 
physics regimes requiring different mindsets. Neverthe- 
less, the variational approach is a useful tool for obtaining 
quantitative insights into the energetics of various phases, 
since it applies directly to the spin- 1/2 problem at hand 
and goes beyond approximate treatments like large-5 and 
mean field. Encouraged by our success for the isotropic 
problem in the field, we apply this tool to the anisotropic 
case, while being critical of the limitations of the varia- 
tional approach. 

We consider the anisotropic triangular lattice antifcr- 
romagnet with J 1 / J < 1, where J and J' are the cou- 
pling constants of horizontal and oblique nearest neigh- 
bour links (see Fig. [1} . The spatially anisotropic wave- 
functions used in this section contain appropriate modi- 
fications to the wavefunctions in Sec. [TXJ 

Permanent constructions: For the wavefunctions from 
Sees. Ill Al III B[ and HID I the localized orbitals used in the 
uud, Y, and V^, C rm now include an additional parameter 
a': 

( 1, r = A 5 
Aoc/ \ ) ~ a ; r = horizontal n.n. of A,-; 
J ' ^ > | —a', r = oblique n.n. of A, ; ^ 
I. 0, otherwise. 

Spiral: Our treatment of the spiral requires separate 
discussion. The 120° spiral generalizes to an incommen- 
surate spiral with wavevector Q&l^ However, in a finite 
sample, periodic boundary conditions would bias against 
the incommensurate order. We can mitigate this effect 
by considering appropriate phase twists at the boundaries 
that accommodate such Q. For computations, it is conve- 
nient to perform a gauge transformation that spreads the 
twist uniformly across the sample; the resulting Hamil- 
tonian is then translationally invariant: 

^twisted = - (trr'e^^'bl br> + h.C.) + ff int , (21) 

(rr f ) 
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where e rr < is the displacement vector from r to r'. The 
twisted Hamiltonian is used only for calculating the in- 
commensurate spiral energies while all other trial ener- 
gies are evaluated using the original Hamiltonian with no 
twist. 

Jastrow factors: To accommodate spatial anisotropy, 
we introduce additional parameters into the nearest 
neighbour and long-range pseudo-potentials as follows: 



u{r, r ) 



w, 



if r and r' are horizontal n.n. 
if r and r' are oblique n.n. 

otherwise. 



(22) 



[a 2 (x-x') 2 + (y-y') 2 ]p/ 2 ' 



2-parton: To obtain spatially anisotropic versions of 
the 2-parton wavefunctions from Sec. Ill HI we allow 
the mean field hopping amplitudes in Eq. (fTBf to be 
anisotropic: 



t 



(») 



ifrandr' are horizontal n.n. 
£'W, if r and r' arc oblique n.n. 



(23) 



We consider the same fluxes and possible site-dependent 
potentials as in Sec. Ill HI For example, at density 
n = 1/3 we allow x \/3 pattern in the chemical po- 
tential. At other densities, we can consider other appro- 
priate CDW patterns. 

One virtue of the 2-parton states is that they con- 
nect naturally to the decoupled chains limit. Indeed, for 
t'(") = 0, sign^ 1 )^ 2 )] = sign[#)] < 0, the trial wave- 
function Eq. (fT5|) on each chain reduces to 



clia 



,XM) 



i-rr(xi-\ hxjw) 



II' 

i<j 



Tr(xi 



(24) 



with p = pi + p2- The first factor gives correct Marshall 
sign for the ID boson problem with hopping < 0. 
(To be more precise, we assume that the chain length 
L is even and choose periodic or antiperiodic boundary 
conditions for the partons depending whether the num- 
ber of bosons M is odd or even.) This is an accurate trial 
state in the full range of boson densities: For n = 1/2, 
with p = 2 it reduces to the ground state of the Haldanc- 
Shastry chain and is a good approximation to the ground 
state of the Heisenberg chain; for n — !> 0, with p = 1 it 
reproduces the nearly free fermion picture of the dilute 
gas of hard-core bosons; for varying n, by adjusting p 
this state can capture varying Luttinger liquid exponents. 
The 2-parton construction can thus provide a starting 
point for exploring what happens when the chains are 
coupled together. As discussed in Sec. Ill HI the parton 
hopping between the chains with vector potentials satis- 
fying Eq. (|19p can roughly capture the interchain boson 
hopping energy, while site-dependent chemical potentials 
can produce candidate CDW states. We will present this 
in some detail for n = 1/3 and n = 1/6. 

ED calculations: To conclude the discussion of our 
anisotropic setups, we describe the supplementary ED 




0.4 0.6 

8 = ( J-J' ) /J 



FIG. 6: Comparison of spiral, permanent-type uud, and 
2-parton trial energies (per bond) for 12 bosons on the 
anisotropic 36-site cluster. The wavefunctions are anisotropic 
generalizations of those constructed in Sec. [TTJ For S < 0.4, 
the optimal 2-parton wavefunction has a non-zero chemical 
potential on sublattice A and provides an alternative real- 
ization of the uud state. For 8 > 0.4, the chemical potential 
optimizes to zero, probably due to large finite-size gap for such 
anisotropy. On a larger 24 x 24 cluster, the chemical potential 
remains non-zero up to S ~ 0.7, leading us to conjecture that 
in the thermodynamic limit the uud phase persists all the way 
to J'/J -> 0. 



calculations on the 36-site cluster with J and J' cou- 
plings. We compute a few lowest eigenvalues in each 
symmetry sector of the Hamiltonian with no twist and 
also eigenvalues in the zero momentum sector of the 
twisted Hamiltonian Eq. (|21[) with varying Q. At a given 
anisotropy and boson density, the minimum of these ED 
energies is taken to be the ground state energy. Our ED 
calculations are restricted to Nb < 12. The variational 
calculations are performed for the same 36-site cluster 
and also for larger systems. 

We now turn to the results of our anisotropic study. 
For illustration, we present two boson densities. 



C. n = i/3 

Figure [5] shows the trial energies of uud, incommen- 
surate spiral, and 2-parton wavefunctions at density 1/3 
on the 36-site cluster. From the isotropic study, it is 
not surprising that the permanent-type uud wavefunc- 
tion remains a good candidate at low anisotropy. As 
mentioned in Sec. Ill HI the 2-parton U1B wavefunction 
constructed with 0/n fluxes through triangles and a lo- 
calizing chemical potential on one sublattice provides an 
alternative realization of the uud state. For S > 0.3, 
the 2-parton energy becomes lower than the permanent- 
type wavefunction but the y/3 x y/3 chemical potential 
(and therefore the uud phase) persists up to 5 0.4. 
Since our permanent wavefunction uses localized orbitals 
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that only extend to nearest neighbour sites, it fails to 
capture longer range correlations in the chain direction 
expected in a more anisotropic system. We would need 
to use more extended orbitals in the permanent, but we 
have not pursued this. On the other hand, the 2-parton 
realization readily accommodates the lattice anisotropy 
via the parton hoppings, Eq. (|23[) . and provides a sim- 
ple way to continue our study of the uud state to larger 
anisotropy. 

In the highly anisotropic region, we obtain good trial 
energies for the 6x6 system using the 2-parton wavefunc- 
tion without the chemical potential. However, if we con- 
sider the low energy cutoff due to the finite cluster size, 
it is clear that the study cannot resolve the true phase in 
the thermodynamic limit. Specifically, in the decoupled 
chains limit, we obtain a very accurate wavefunction for 
two bosons on a 6-sitc chain by using antiperiodic bound- 
ary conditions for the partons, cf. Eq. (f2"4"]l . The corre- 
sponding parton spectrum nicely accommodates two par- 
ticles and has a large finite-size gap to next levels, which 
persists up to moderate interchain couplings. While our 
2-parton state by virtue of good fluxes naturally builds in 
good interchain exchange correlations, we cannot resolve 
the thermodynamic phase (e.g., the development of the 
v/3 x V3 CDW) if the relevant energy scale is much lower 
than the finite-size gap. 



Spiral 




V ■» 




2p (u>0) ■ 
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0.2 0.4 0.6 0.8 1 

5 = ( J-J' ) /J 

FIG. 7: Comparison of spiral, commensurate VByi_j a strow, 
and 2-parton trial energies as a function of anisotropy for 6 
bosons on the 6x6 cluster. The Ei v values are obtained 
with higher chemical potential on one sublattice for S < 0.3 
(in agreement with this state trying to capture 3-sublattice 
features near small anisotropy in this small sample). We 
also performed a much larger study at n = 1/6 comparing 
commensurate and incommensurate V states, incommensu- 
rate spiral, and 2-parton states; from this study, the region 
of commensurate V is actually quite small, while the incom- 
mensurate V dominates over the spiral over a range 5 < 0.3 
(see text for more details) . 



To determine how far the uud phase might extend into 
the anisotropic region, we repeat the 2-parton calculation 
on a large 24 x 24 cluster and find that the s/3 x -y/3 
chemical potential surprisingly remains non-zero up to 
S ps 0.7. We also check that the parton spectrum for the 
optimal parameters is fully gapped and is connected to 
the strongly gapped CDW limit, so the trial wavefunction 
is indeed a valid charge-ordered Mott insulator of bosons 
as discussed in Sec. Ill HI We thus conclude that the 
uud state persists to rather strong anisotropy, albeit the 
CDW order becomes progressively weaker. Interestingly, 
Ref. [13 would predict the same x CDW order m 
the nearly decoupled chains limit at density n = 1/3. 
Combining with our variational work, this suggests that 
the uud phase may in fact extend continuously from 5 = 
up to 5 = 1 (See Fig. [TU]). A rigorous confrontation 
to this conjecture could be provided, for example, by 
a systematic DMRG study of 3 x L ladders at density 
n = 1/3 varying J'/ J from 1 to and monitoring the 
evolution of the v3 x vo charge order. 

While our study agrees with the observed plateau in 
Cs2CuBr4 (S w 0.3), it contradicts the absence of the 
plateau in CS2CUCI4 (5 « 0.66). It is likely that residual 
interactions (e.g. such as Dzyaloshinskii-Moriya) have to 
be added to the Heisenberg model in order to describe 
the latter material ^ and they can change the energetics 
balance against the (very weak) uud state in this highly 
anisotropic system. 



D. n = l/6 

Figure [7] shows the trial energies of commensurate 
Veiji-jastrow, incommensurate spiral [treated as described 
around Eq. (|21[) ]. and 2-parton wavefunctions for 6 
bosons on the 6x6 cluster. At low anisotropy, the 
Veiji-jastrow state is a good candidate. For increas- 
ing anisotropy, a change to the incommensurate spiral 
is observed, which eventually loses to the "quasi- Id" 
phase represented by the 2-parton wavefunction with zero 
chemical potential. 

In the highly anisotropic region, the figure shows re- 
markable agreement between the ED and the 2-parton 
energies, where we impose uniform 7r/6 and 57r/6 flux per 
triangle for the d\ and di partons respectively (Chern- 
Simons state described in Sec. Ill Hp . Despite the ex- 
cellent agreement, we simply conclude that the highly 
anisotropic region is strongly dominated by quasi- Id 
physics and finite size effects. Specifically, in the ED 
calculation on the 36-sitc cluster with 6 bosons, we find 
a non-degenerate ground state and a relatively large ex- 
citation energy gap. We interpret this as follows. In the 
decoupled chains limit, each chain contains one boson; 
for such a segment of length L = 6, one expects a non- 
degenerate ground state with a large excitation gap due 
to finite size. This gap persists as the chains are coupled, 
particularly because of some frustration present in the 
triangular lattice geometry. 

We can similarly rationalize all our ED observations at 
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other densities in the highly anisotropic limit. For exam- 
ple, for 7 bosons on the 36-site cluster, one of the chains 
now contains two bosons, and the ground state of the 
decoupled chains Hamiltonian is 6-fold degenerate due 
to 6 possible ways of choosing this chain. The finite-size 
gaps "protect" this situation until the interchain coupling 
J' becomes sufficiently large. Such observations on the 
ED spectra show serious limitations of the small system 
study in the anisotropic model. Going over all ED data 
for N b < 12, we conclude that 5 > 0.5 regime can be ra- 
tionalized as such weakly-coupled finite chains, with no 
clear resolution of the ultimate state. This is labeled as 
"quasi-ld" region in Fig. [9j 

One of the goals of the 6x6 study was to have ED 
reference for our trial states. Having achieved some con- 
fidence in the good energetics of these states (despite 
their limitations) , we now want to discuss variational re- 
sults for larger sizes. Specifically, on the 36-site cluster, 
we have not considered the possibility of incommensu- 
rate V state: while we know how to accommodate the 
incommensurate spiral state, we do not have similar con- 
struction for the incommensurate coplanar state. On the 
36-site cluster, we see that incommensuration becomes 
important for 5 > 0.2. In fact, as we discuss below, we 
think that for this density the V state is probably in- 
commensurate already for smaller anisotropy, but also 
extends to larger anisotropy in the competition against 
the spiral. 



E. Incommensurate V versus spiral study at low to 
intermediate boson densities 

In this section, we focus on the high-field regime where 
the incommensurate V and spiral are the main competing 
candidates. First, we briefly describe the relevant physi- 
cal picture. Beginning with the near-saturation limit, we 
consider a gas of free bosons hopping on the triangular 
lattice with the following kinetic energy spectrum: 

e k = J cos{k x ) + 2 J' cos (^fj cos ■ (25) 

The band minima occur at Q = ± (Q x , 0) with 

Q x = 2 acos(-J'/2J). (26) 

A condensation of bosons at these points gives rise to a 
degenerate manifold of states spanned by 

{(^H&l^-^IO) ; m = 0,l,...,N b }. (27) 

At low densities, the degeneracy is lifted by nearest neigh- 
bour repulsion. To see how this happens, we expand the 
interaction in terms of the two dominant spectral modes 



and then replace the operators by c-number a 11 ' 12 : 

b r ~ e i<3? b$ + <T®* (28) 

H in t = Jrr ' b r b r b t' br ' ( 29 ) 

(rr') 

~ (j + 2j')(|6qI 2 + I&_qI 2 ) 2 

+ 2v \b \ 2 \b^\\ (30) 
v = Jcos{2Q x ) + 2J'cos(Q x ). (31) 

The effect of nearest neighbour repulsion is determined 
by the sign of v. For J' / J < 0.39, v is positive and the 
incommensurate spiral (e.g., I&^l ^ and \b_^\ = 0) 
wins. For 0.39 < J'/ J < 1.59, v is negative and the 
incommensurate V ordering = |&_q|) becomes more 
stable. The prediction from this approximate treatment 
is consistent with the spiral phase found in the nearly 
decoupled chains limit near saturation and in the highly 
anisotropic dilute boson study for the Cs2CuCi4 j 10 ' 12 this 
is also consistent with the coplanar phase found in the 
isotropic dilute boson studyjii 

In the above discussion, we have neglected the effect of 
hard core interaction. Intuitively, this should be more im- 
portant at higher density: The role of the hard core con- 
straint is to prevent two bosons already in nearest neigh- 
bour contact from further occupying the same site, while 
at low density such contacts are avoided due to nearest 
neighbor repulsion. To see whether the hard core inter- 
action favors the V or spiral phase, we expand the on-site 
repulsion energy in terms of the two spectral modes: 

{blKf ~ (|fo Q -| 2 + |6_ Q H 2 ) 2 + 2 |fe Q H 2 Kq\ 2 - (32) 

From the positive sign in the second term, which dislikes 
the V, we may expect the boundary between the V and 
spiral phases to shift in favor of the spiral phase as the 
density increases. 

To address the competition between these two phases 
quantitatively at finite density, we implement a varia- 
tional study on m x n rectangular clusters such that 
a fitting wavevector Q x = 2np/m is close to the spiral 
wavevector at each anisotropy (with appropriate integers 
p). For the spiral phase, we use the same anisotropic 
wavefunction described earlier. A candidate wavefunc- 
tion for the incommensurate V phase is constructed as 
follows: 

llfcO = (^ + e~ ia b^y > \0) 

= ^cos(Q-f+a) b{\ |0>. (33) 

Note that for incommensurate Q the relative phase be- 
tween b - and b^ - is not fixed, which we indicated with 
Q -Q 

a. This is not important in an infinite system since Q ■ r 
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FIG. 8: Variational phase boundary between V and spiral (S) 
phases obtained from calculations on cluster sizes 18 x 18, 51 x 
20, 27x20, 44x20, 47x20, 38x10, 39x10, 53x10, 60x10, and 
99 x 10, where the corresponding anisotropy increases from 
to 0.9 in equal intervals; the sizes are chosen so as to best 
accommodate the classical wavevector, Eq. (|26|l . The V phase 
remains commensurate at 5 = 0.1 for n > 1/5. Note that at 
high anisotropy and particularly with increasing density, we 
find that the spiral loses to the 2-parton wavefunction, which 
we interpret as quasi- ID dominated regime. 



visits all phases. On the other hand, for commensurate 
Q = (47r/3, 0), a = and tt/2 correspond to distinct V 
and ^P-type phases discussed in the isotropic case; both 
can be viewed as "parent" states for the incommensu- 
rate coplanar phase, but we will continue referring to the 
latter as V-typc. 

For ease of implementation, we use the same trans- 
lationally invariant pseudopotentials given in Eq. (|22p. 
The V state has an incommensurate density wave and in 
principle allows more complicated pseudopotentials, so 
this choice probably biases slightly in favor of the spiral 
which has uniform boson density. In all other respects, 
the physical setting and the variational freedom are very 
similar in our realizations of the spiral and V states, and 
we think this study provides a fair comparison between 
the two phases even if the Q x may be slightly off and the 
Jastrow pseudopotentials are not the most general. 

Fig.[5]shows the result of our incommensurate V versus 
spiral variational study. The boundary between the two 
phases qualitatively agrees with our earlier argument, 
suggesting that the hard core repulsion is comparatively 
less important at low density. We note that the obtained 
trial energies of the V and spiral states are quite close 
(particularly at low density) , hence the exact location of 
the phase boundary should not be taken as definitive. 
Furthermore, the simple pscudopotcntial is clearly not 
optimal in the highly anisotropic regime, and eventually 
our spiral loses to the 2-parton states. A more rigorous 
V versus spiral variational study can be pursued by in- 
troducing more variational parameters into the Jastrow 
factor and employing systematic wavefunction optimiza- 
tion methods^ The variational result that the coplanar 
phase extends to large anisotropy for fields near satu- 
ration is also in agreement with a recent dilute boson 
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FIG. 9: Variational phase diagram for the anisotropic 6x6 
cluster in magnetic field. In this study, we exclude incom- 
mensurate versions of V and Y states. Due to finite size limi- 
tation, certain regions are marked as unresolved. In addition, 
the region S > 0.5 is strongly dominated by quasi-ld physics, 
particularly for this small cluster study (see discussions in 
text). 

calculation extended to all J'/ J^i 

From the present results, we make an interesting ob- 
servation that at the anisotropy relevant for Cs2CuBr4, 
5 w 0.3, the transition occurs at density somewhere be- 
tween n = 1/5 (magnetization 0.6 of saturation) and 
n = 1/6 (magnetization 2/3 of saturation). The V phase 
occupies the region near saturation, while the spiral oc- 
curs at lower magnetizations. Thus, if the Hciscnberg 
model is an adequate description, some of the features in 
the high field phase diagram of Cs2CuBr4 may be due to 
the competing umbrella-type and coplanar states^ 



F. Summary of anisotropic study 

Figure [5] summarizes a variational phase diagram ob- 
tained for the 36-site cluster considering all boson densi- 
ties. We label certain parts of the diagram with question 
marks or broken lines to indicate these regions as unre- 
solved or less reliable. The figure shows the uud phase 
extending relatively far into the anisotropic region. On 
both sides of the uud phase, the commensurate copla- 
nar phases remain stable over the incommensurate spi- 
ral for certain ranges of the anisotropy. As the spatial 
anisotropy biases against the commensurate states, the 
actual V and Y regions are expected to be wider if the 
wavefunctions are generalized to incommensurate ver- 
sions. However, we exclude such extensions since they 
could not be accommodated on the 6x6 cluster. For 
5 > 0.5, our 2-parton trial energies are generally very 
good. However, we think that this only indicates the 
onset of quasi-ld physics and strong finite-size effects as 
discussed earlier for the specific densities. 
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FIG. 10: Schematic phase diagram for the anisotropic trian- 
gular lattice in magnetic field, combining 6x6 study as well as 
larger cluster studies. We suggest the possibility that the uud 
plateau extends across the entire range of anisotropy. The V 
phase is commensurate (see Fig. [2} near the isotropic axis 
and the plateau but becomes incommensurate at moderate 
anisotropy and higher fields. The highly anisotropic region 
is not well resolved in our variational study. The low field 
regime is also not studied thoroughly in this work, while pre- 
vious studie d 20 ' 23 ' 25 suggest spin liquid state along h = and 
8 > 0.2. 



Figure [10] shows a schematic phase diagram based on 
the 6x6 anisotropic study as well as studies on larger 
clusters. Here, we address a number of unresolved re- 
gions in Fig. 121 limits of the uud plateau, the boundary 
between V and spiral, and the boundary between com- 
mensurate and incommensurate V. We find that the uud 
phase extends much futher and may be even to all S. 
Also, a significant portion of the phase diagram at high 
fields is occupied by the incommensurate V phase. 

We note that our work does not rule out other in- 
commensurate phases found in the recent study by Al- 
icea. et ai£ For example, we were not able to come up 
with a good implementation of the incommensurate ex- 
tension of the Y state. On the other hand, we did try Bijl- 
Jastrow-type wavefunctions for distorted umbrella states 
discussed in Ref. 0, which are commensurate supersolids 
with incommensurate spiral phase angles. On the 36-site 
cluster, these trial states optimized to the incommensu- 
rate spiral with uniform boson density, but we have not 
explored this thoroughly on larger clusters. Overall, our 
results are more conclusive at low densities and much less 
at high densities between 1/3 and 1/2. 



IV. SUMMARY AND DISCUSSION 

We studied the Heisenberg antifcrromagnct on the 
spatially anisotropic triangular lattice in the field from 
a variational perspective. On the isotropic lattice, we 



constructed a very simple and physically transparent 
permanent-type wavefunction for the uud state at density 
1/3. This is a Mott insulator of bosons where we accu- 
rately included small charge fluctuations by using appro- 
priate localized boson orbitals. The remarkable trial en- 
ergy suggests that such approach may be useful in other 
Mott insulator contexts. Next, we obtained natural ex- 
tensions to nearby V and Y supersolid phases respec- 
tively for n < 1/3 and n > 1/3, where the physics re- 
mains strongly influence by the proximity to n = 1/3. 
By connecting to a Bijl-Jastrow-type candidate wave- 
function at low density, the coplanar V phase extends 
to all n < 1/3 (i.e., up to the saturation field in the spin 
model language). Note, however, that at very low density 
another coplanar state (vP-type) is expected to be very 
closed and we cannot resolve between the two. On the 
higher density side of the plateau (i.e., at lower fields), 
the permanent-type Y wavefunction performs well near 
the plateau but narrowly loses to the Huse-Elser spiral 
candidate at densities close to half- filling (zero field). The 
latter result is consistent with other recent worksi 21 ' 24 ' 26 
The success of our isotropic study encouraged us to 
extend it to the anisotropic lattice. At density n = 1/3, 
we begin with the permanent-type realization of the uud 
and then connect to a conceptually similar but techni- 
cally different 2-parton realization at higher anisotropy. 
Surprisingly, we found that the uud phase remains stable 
over a large range of anisotropy. In conjunction with the 
same CDW phase found in the decoupled chains limitj 1 ^ 
we suggest that the uud phase may in fact extend across 
the entire range of anisotropy. This conjecture can be 
tested more rigorously using a DMRG study on finite- 
width strips. 

In the low boson density region (i.e., at high fields), 
the Bijl-Jastrow-type V commensurate supersolid wave- 
function is smoothly connected to the incommensurate 
version. This state competes with the incommensurate 
spiral, and we can accurately compare the two. We found 
that the incommensurate V state has lower energy in a 
large region of the phase diagram, extending up to a fairly 
large value of anisotropy in the very dilute regime (i.e., 
close to the saturation fields). On the other hand, the V 
phase remains commensurate near the isotropic axis and 
the plateau. 

In the high density regime (i.e., at low fields), we at- 
tempted to construct an incommensurate Y candidate 
using a Bijl-Jastrow-type wavefunction but found that 
this construction performs poorly. This low field region 
at moderate to high lattice anisotropy calls for more com- 
prehensive investigation. 

One of the goals we had was to explore possible new 
plateaus in the high field regime of Cs2CuBr4. We have 
learned that the phase diagram is already very rich even 
without considering any additional plateaus. Neverthe- 
less, for several densities such as 1/6, 2/9, and 1/4 we 
implemented permanent-type wavefunctions for various 
proposed CDW from Ref. Qj as well as for some addi- 
tional stripe-like ordcrings, and inevitably found that ei- 
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ther V or spiral has lower energy. Our earlier uud study 
showed that the 2-parton construction can also be useful 
for studying CDW phases; however, similar implemen- 
tations at the above densities again failed to reveal any 
stable charge ordering. This suggests that any such or- 
der, if present at all, is likely to be very weak. 

One of the findings from our study is that for the 
Cs2CuBr4 anisotropy, the system is in the incommensu- 
rate coplanar phase close to the saturation fields^ and 
there may be a transition to the non-coplanar spiral state 
at lower fields; this could be responsible for one of the fea- 
tures in the Cs2CuBr4 experiment. We cannot exclude 
other more complex cascade of phases. Furthermore, ad- 
ditional residual interactions not treated here may be im- 
portant for understanding the phases of Cs2CuBr4 in the 
field. This remains a fascinating open problem. 



Appendix A: Motivation for V pelm wavefunction for 
n < 1/3, Eq. (TT^I 

We begin with the uud state with N/3 bosons, Eq. ((6}, 
and put Nh = N/3 - N b holes in a "hole orbital" <f>h(R), 




For a boson configuration 

\v) = 4 1 ---CJ0) , (A2) 
we obtain an amplitude, 



/ 0i oc M 



J N/3 
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D N/3 



Aloc 3 (-Ri) 
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\<f>f c (R Nh ) ... 4>%(RN h )J 



4>? c (r Nb ) 

Cl 



CN/3 



CN/3 J 



(A3) 



(tj|¥) = Y! MRi) ■ --MRnJ Perm 

Ri...R Nh 



The primed sum indicates that R\, . . . 1 Rn h need to be 
different from each other and from all n , . . . , rjv 6 • Close 
to the 1/3 plateau, the density of holes is small, and 
we can approximately replace the restricted sum by an 
unrestricted sum. Performing independent summations 
over R\, . . . , Rn h gives the last expression in the form of 
a single permanent, where 

Cj = ymr)^ c ( r ) ( A4 ) 

R 

is an "overlap" of the 4>h and 4> x ° c orbitals. On physics 
grounds, the hole orbital 4>h needs to respect the symme- 
tries of the uud state (cfih = const over the A sublattice). 
In this case, Cj is independent of j, and up to a nor- 
malization constant we can replace all matrix elements 
in the last Nh rows by 1. The approximate single per- 
manent form is a valid variational wavefunction by itself, 
and this is the state we use in the main text and call 

Appendix B: Correlation functions of 
permanent-type states 

We calculate the order parameters and correlation 
functions in the permanent-type trial states, since we 
do not have much experience with such wavefunctions 



I 

for Mott insulators or supcrsolids. The non-permanent 

V and spiral trial states are more obvious constructions, 
and therefore omitted. 

For Y, uud, and V'pcrm trial states, a three sublat- 
tice modulation is observed in the CDW order param- 
eter (n r ). As expected, the density structure factor 
(n- q n q ) reveals sharp peaks near the reciprocal vectors 
Q = ±(4p0) for all three trial states. 

Figure [TT] shows the correlation functions (b\b r i) 
between two sites for these trial states. The uud state 
has rapidly decaying correlations between any two sites, 
which is expected in this Mott insulator state. For the 

V supersolid state, the correlations decay rapidly when 
at least one site lies on the higher density sublattice 
A (i.e., as if this sublattice is Mott- insulating), while 
they are long-ranged when both sites reside on the BC 
honeycomb sublattice. The signs of the correlations are 
positive for all pairs of B-B or C-C sites, and negative 
for all B-C pairs, which is consistent with the Y spin 
order shown in Fig. [2j Finally, for the V vcrm trial state, 
long-ranged correlation exists between any two sites on 
the lattice. The signs are negative for all pairs of A-B 
or A-C sites, and positive for all B-C pairs (as well as 
A- A, B-B, and C-C pairs), which is consistent with the 

V spin order shown in Fig. [2] Thus, we have verified our 
intuition about the physical properties of these states. 
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Y: n=0.38 



uud: n=l/3 



Vperm : n=0 . 2 9 




FIG. 11: Boson correlation (6j6 r /) as a function of real space distance \r — r'\. Top row: r belongs to sublattice A (where 
tia > riB = nc) while r' includes sites on all three sublattices. Bottom row: r belongs to sublattice B while r' includes sites 
on sublattices B and C. Hexagonal clusters with 84, 144, and 48 sites are respectively used for Y (left column), uud (middle), 
and Vpcrm (right) calculations. 
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